--- title: OCO2 - Reproduce F. Chevallier's figures keywords: fastai sidebar: home_sidebar ---
!apt-get install libgeos-3.5.0
!apt-get install libgeos-dev
!pip install https://github.com/matplotlib/basemap/archive/master.zip
from google.colab import drive
drive.mount('/content/drive')
Using Data from OCO-2 Satellite, issued by the NASA.
We here try to reproduce the results from this paper by F. Chevallier, trying to "[observe] carbon dioxide emissions over China's cities with the Orbiting Carbon Observatory-2".
//TODO: Explanation
import pandas as pd
import numpy as np
from numpy import exp, loadtxt, pi, sqrt
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
from mpl_toolkits.basemap import Basemap #Imported directly from the github repository
Parameters:
def draw_map(data, x="longitude", y="latitude", c="xco2", lon_min=-180, lon_max=180, lat_min=-90, lat_max=90, size_point=1, frontier=False):
plt.figure(figsize=(15, 10), edgecolor='w')
m = Basemap(llcrnrlat=lat_min, urcrnrlat=lat_max, llcrnrlon=lon_min, urcrnrlon=lon_max)
m.shadedrelief()
parallels = np.arange(-80.,81,10.)
m.drawparallels(parallels,labels=[False,True,True,False])
meridians = np.arange(10.,351.,20.)
m.drawmeridians(meridians,labels=[True,False,False,True])
normal = matplotlib.colors.LogNorm(vmin=data[c].min(), vmax=data[c].max())
m.scatter(data[x], data[y], c=data[c], cmap=plt.cm.jet, s=size_point, norm=normal)
if (frontier):
m.drawcountries(linewidth=0.5)
m.drawcoastlines(linewidth=0.7)
plt.show()
Parameters:
Return:
from math import radians, cos, sin, asin, sqrt
def haversine_formula(lon1, lat1, lon2, lat2):
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
r = 6371 # Radius of earth in kilometers. Use 3956 for miles
return c * r
def haversine(row, lat, lon):
return haversine_formula(lon, lat, row['longitude'], row['latitude'])
Sample data can be accessed freely on the NASA Database, among other open data from several NASA sattelites.
We will be using a CSV of peaks detected by Benoit Courty here.
data = pd.read_csv("https://courty.fr/OCO2/result_doco2_1808-o22158.csv", sep=",", index_col=0)
data.head()
data.describe()
We spot some negative sigma and amplitude parameters.
data['sigma'] = data['sigma'].apply(abs)
data['amplitude'] = data['amplitude'].apply(abs)
data.head(15)
data.describe()
To convert the sounding_id into a datetime variable data:
from datetime import datetime
def to_date(a):
return datetime.strptime(str(a), '%Y%m%d%H%M%S%f')
data['date'] = data['sounding_id'].apply(to_date)
data.head()
import scipy.stats as stats
df = data[:5]
for sigma in df['sigma'][:5]:
std = abs(sigma)
x = np.linspace(-15, 15, 100)
y = stats.norm.pdf(x,0,std)
plt.plot(x,y, color='coral')
plt.grid()
plt.show()
from plotly.figure_factory import create_quiver
import plotly.graph_objects as go
import plotly.express as px
fig = px.scatter_mapbox(data, lat="latitude", lon="longitude",
hover_name="date", hover_data=["amplitude", "sigma"],
color_discrete_sequence=["red"], zoom=2, height=700)
fig.update_layout(mapbox_style="open-street-map")
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
fig.show()
fig.write_html("./peaks.html")
import folium
folium_map = folium.Map([43, -100], zoom_start=4,tiles="CartoDB dark_matter")
for index, row in df_res.iterrows():
radius = row["amplitude"]/20
color="#E37222" # tangerine
tooltip = "["+str(round(row['latitude'],2))+" ; "+str(round(row['longitude'],2))+"]"
folium.CircleMarker(location=(row["latitude"],
row["longitude"]),
radius=radius,
color=color,
tooltip=tooltip,
fill=True).add_to(folium_map)
folium_map
!pip install geopandas
import os
import glob
import numpy as np
import pandas as pd
input_dir = 'http://courty.fr/OCO2/'
# http://courty.fr/OCO2/peak_data-f_oco2_1808-o_21733-si_2018080216413233.json
def load_one_peak_data(sounding_id, df_all_peak):
df_param = df_all_peak.query("sounding_id==@sounding_id")
param_index = df_param.index[0]
url = input_dir + "peak_data-f_oco2_1808-o_"+str(df_param.loc[param_index, 'orbit'])+"-si_"+sounding_id+".json"
print(url)
df_peak = pd.read_json(url)
gaussian_param = {
'slope' : df_param.loc[param_index, 'slope'],
'intercept' : df_param.loc[param_index, 'intercept'],
'amplitude' : df_param.loc[param_index, 'amplitude'],
'sigma': df_param.loc[param_index, 'sigma'],
'delta': df_param.loc[param_index, 'delta'],
'R' : df_param.loc[param_index, 'R'],
}
return df_peak, gaussian_param
'''
x : the data input value
m : the slope of the data
b : the intercep of the data
A : Amplitude de la courbe
sig : sigma / écart type de la courbe
'''
def gaussian(x, m, b, A, sig):
return m * x + b + A / (sig * (2 * np.pi)**0.5) * np.exp(-x**2 / (2*sig**2))
def plot_peak(df_peak, gaussian_param):
x = df_peak['distance']
y = df_peak['xco2']
y = y - gaussian_param['slope'] * x - gaussian_param['intercept']
plt.scatter(x, y, c=y, s=3, label='sounding')
plt.plot(x, gaussian(x, m=0, b=0, A=gaussian_param['amplitude'], sig=gaussian_param['sigma']), 'r', label='fit')
plt.legend()
plt.title('OCO 2 data')
plt.xlabel('Distance')
plt.ylabel('CO²')
plt.show()
df_res = pd.read_csv(input_dir + "result_for_oco2_1808.csv", sep=",")
df_peak, gaussian_param = load_one_peak_data("2018083122324108", df_res)
print(gaussian_param)
plot_peak(df_peak, gaussian_param)
import folium
folium_map = folium.Map([43, -100], zoom_start=4,tiles="CartoDB dark_matter")
for index, row in df_res.iterrows():
radius = row["amplitude"]/20
color="#E37222" # tangerine
tooltip = "["+str(round(row['latitude'],2))+" ; "+str(round(row['longitude'],2))+"]"
# popup = folium.Popup(max_width=450).add_child(folium.Vega(peak, width=450, height=250))
# url = input_dir + "peak_data-f_oco2_1808-o_"+str(int(row['orbit']))+"-si_"+str(int(row['sounding_id']))+".json"
# print(url)
# peak = pd.read_json(url)
folium.CircleMarker(location=(row["latitude"],
row["longitude"]),
radius=radius,
color=color,
tooltip=tooltip,
# popup=popup,
fill=True).add_to(folium_map)
folium_map
# Create quiver figure
fig = create_quiver(data.longitude, data.latitude, data.windspeed_u, data.windspeed_v,
scale=.25,
arrow_scale=.4,
name='quiver',
line_width=1)
fig.add_trace(go.Scatter(x=data.longitude, y=data.latitude,
mode='markers',
marker_size=2,
name='points'))
fig.update_layout(mapbox_style="open-street-map")
fig.update_layout(margin={"r":0,"t":0,"l":0,"b":0})
fig.show()
fig.write_html("./winds.html")
path_Edgar_2018 = "/content/drive/My Drive/Data For Good - S7: Ministère & O-CO2/data/edgar/CO2_emissions_Edgar_2018_v3.csv"
edgar = pd.read_csv(path_Edgar_2018, sep=",", index_col=0)
edgar.head()
# We only keep the 75% highest registered emitters
edgar_top = edgar[edgar['CO2 classification encoded'] > 2]
edgar_top.describe()
edgar_top['CO2 classification encoded'][round(0.95*dataframe_list_v1.shape[0]):] = 'Very high'
folium_map = folium.Map([43, -100], zoom_start=4,tiles="CartoDB dark_matter")
for index, row in df_res.iterrows():
radius = row["amplitude"]/20
color="#E37222" # tangerine
tooltip = "["+str(round(row['latitude'],2))+" ; "+str(round(row['longitude'],2))+"]"
# popup = folium.Popup(max_width=450).add_child(folium.Vega(peak, width=450, height=250))
# url = input_dir + "peak_data-f_oco2_1808-o_"+str(int(row['orbit']))+"-si_"+str(int(row['sounding_id']))+".json"
# print(url)
# peak = pd.read_json(url)
folium.CircleMarker(location=(row["latitude"],
row["longitude"]),
radius=radius,
color=color,
tooltip=tooltip,
# popup=popup,
fill=True).add_to(folium_map)
for index, row in edgar_top.iterrows():
radius = 1
color="#3186CC" # blue
tooltip = "["+str(round(row['latitude'],2))+" ; "+str(round(row['longitude'],2))+"]"
folium.CircleMarker(location=(row["latitude"],
row["longitude"]),
radius=radius,
color=color,
tooltip=tooltip,
fill=True).add_to(folium_map)
folium_map
path_plants = "https://raw.githubusercontent.com/dataforgoodfr/batch7_satellite_ges/master/dataset/CO2_emissions_centrale.csv"
plants = pd.read_csv(path_plants, sep=",", index_col=0)
plants.head()
plants.describe()
peaks_and_plants_dark = folium.Map([43, 0], zoom_start=4,tiles="CartoDB dark_matter")
for index, row in plants.iterrows():
radius = 0.1
color="#3186CC" # blue
tooltip = "["+str(round(row['latitude'],2))+" ; "+str(round(row['longitude'],2))+"]"
emit = str(round(row['estimated_generation_gwh'],2))
popup_html="""<h4>"""+tooltip+"""</h4>"""+row['country_long']+"""<p>Emission 2018 (est): """+emit+""" GWh</p>"""
popup=folium.Popup(popup_html, max_width=450)
folium.CircleMarker(location=(row["latitude"],
row["longitude"]),
radius=radius,
color=color,
tooltip=tooltip,
popup=popup,
fill=True).add_to(peaks_and_plants_dark)
for index, row in data.iterrows():
radius = row["amplitude"]/20
color="#E37222" # tangerine
tooltip = "["+str(round(row['latitude'],2))+" ; "+str(round(row['longitude'],2))+"]"
sounding = str(row['sounding_id'])
date = str(row['date'])
orbit = str(row['orbit'])
popup_html="""<h4>"""+tooltip+"""</h4>"""+date+"""<p>sounding_id: """+sounding+"""</br>orbit: """+orbit+"""</p>"""
popup=folium.Popup(popup_html, max_width=450)
folium.CircleMarker(location=(row["latitude"],
row["longitude"]),
radius=radius,
color=color,
tooltip=tooltip,
popup=popup,
fill=True).add_to(peaks_and_plants_dark)
peaks_and_plants_dark